rm(list=ls())
data <-read.csv("Data.csv")

#install.packages("doBy")
#install.packages("gmodels")
#install.packages("ggplot2")
#install.packages("xtable")
#install.packages("XLConnect")
#install.packages("xlsx")
#install.packages("aod")
#install.packages("MASS")
#install.packages("mvtnorm")
#install.packages("xtable")
#install.packages("stargazer")
#install.packages("dplyr")
#install.packages("descr")
#install.packages("DescTools")
#install.packages("ggpubr")
#install.packages("purrr")
#install.packages("forcats")
#install.packages("showtext")

#install.packages("extrafont")
#library(extrafont)
#font_import(prompt = FALSE)
#loadfonts()

library(doBy) 
library(plyr)
library(dplyr)
library(gmodels) 
library(ggplot2)
library(xtable)
library(XLConnect)
library(xlsx)
library(sandwich)
library(lmtest)
library(aod)
library(MASS)
library(xtable)
library(stargazer)
library(descr)
library(DescTools)
library(ggpubr)
library(scales)
library(rlang)
library(forcats)
library(showtext)


font_add_google("EB Garamond")
windows()

showtext_auto()



# Main Findings #

## Figure 2

get_mean_sd <- function(vector){
  mean <- mean(vector, na.rm = T)
  sd <- sd(vector, na.rm = T)
  sem <- sd/(sqrt(length(vector)))
  output <- paste(round(mean, 3), " (", round(sem, 3), ")", collapse = "", sep = "")
  return(output)
}

PerceivedSafety_table <- group_by(data) %>%
  summarize(Base_Safe = get_mean_sd(Base_Safe),
            Brigade_Safe = get_mean_sd(Brigade_Safe),
            Company_Safe = get_mean_sd(Company_Safe),
            Stockpile_Safe = get_mean_sd(Stockpile_Safe),
            Destroyers_Safe = get_mean_sd(Destroyers_Safe),
            Flybys_Safe = get_mean_sd(Flybys_Safe))

PerceivedSafety_table

SafetyData <- data.frame(
  name=c("Fighting Force: Permanent U.S. Base", "Fighting Force: U.S. Army Brigade","Offshore: Stockpiling Supplies", "Offshore: U.S Navy Destroyer", "Transient Demonstration: Bomber/Fighter Patrols", "Tripwire: 240 U.S. Soldiers"),
  mean=c(4.217, 4.093, 3.364, 3.209, 3.023, 2.93),
  se=c(.098,.087,.127,.123,.119,.136)) 

tiff("Figure2_Safety.tif", res=100, height=10, width=10, units="in")
SafetyData %>%
  arrange(mean) %>%
  mutate(name = factor(name, levels=c("Fighting Force: Permanent U.S. Base", "Fighting Force: U.S. Army Brigade", "Offshore: Stockpiling Supplies", "Offshore: U.S Navy Destroyer", "Transient Demonstration: Bomber/Fighter Patrols", "Tripwire: 240 U.S. Soldiers"))) %>%
  ggplot(aes(x = name, y = mean)) + geom_point() +
    geom_errorbar(aes(ymin = mean - 1.96 *se, ymax = mean + 1.96 * se), width = .2) +
    theme_minimal(base_size=12, base_family="EB Garamond") + xlab("Reassurance Measure") + ylab("Mean Perceived Safety") +
    theme(axis.text.x = element_text(hjust = 1)) + coord_flip()
dev.off()


## Figure 3

PerceivedLikelihood_table <- group_by(data) %>%
  summarize(Base_Resolve = get_mean_sd(Base_Resolve),
            Brigade_Resolve = get_mean_sd(Brigade_Resolve),
            Company_Resolve = get_mean_sd(Company_Resolve),
            Stockpile_Resolve = get_mean_sd(Stockpile_Resolve),
            Destroyers_Resolve = get_mean_sd(Destroyers_Resolve),
            Flybys_Resolve = get_mean_sd(Flybys_Resolve))

PerceivedLikelihood_table

LikelihoodData <- data.frame(
  name=c("Fighting Force: Permanent U.S. Base", "Fighting Force: U.S. Army Brigade", "Tripwire: 240 U.S. Soldiers", "Offshore: Stockpiling Supplies", "Offshore: U.S Navy Destroyer", "Transient Demonstration: Bomber/Fighter Patrols"),
  mean=c(3.978, 3.93, 3.395, 3.311, 3.205, 3.209),
  se=c(0.165, 0.145, 0.129, 0.137, 0.14, 0.149)) 

tiff("Figure3_Willingness.tif", res=100, height=10, width=10, units="in")
LikelihoodData %>%
  arrange(mean) %>%
  mutate(name = factor(name, levels=c("Fighting Force: Permanent U.S. Base", "Fighting Force: U.S. Army Brigade", "Tripwire: 240 U.S. Soldiers", "Offshore: Stockpiling Supplies", "Offshore: U.S Navy Destroyer", "Transient Demonstration: Bomber/Fighter Patrols"))) %>%
  ggplot(aes(x = name, y = mean)) + geom_point() +
  geom_errorbar(aes(ymin = mean - 1.96 *se, ymax = mean + 1.96 * se), width = .2) +
  theme_minimal(base_size=12, base_family="EB Garamond") + xlab("Reassurance Measure") + ylab("Mean Perceived Likelihood that U.S. Will Aid Ally in Defense") +
  theme(axis.text.x = element_text(hjust = 1)) + coord_flip()
dev.off()

## Figure 4

PerceivedCapability_table <- group_by(data) %>%
  summarize(Base_Capability = get_mean_sd(Base_Capability),
            Brigade_Capability = get_mean_sd(Brigade_Capability),
            Company_Capability = get_mean_sd(Company_Capability),
            Stockpile_Capability = get_mean_sd(Stockpile_Capability),
            Destroyers_Capability = get_mean_sd(Destroyers_Capability),
            Flybys_Capability = get_mean_sd(Flybys_Capability))

PerceivedCapability_table

CapabilityData <- data.frame(
  name=c("Fighting Force: Permanent U.S. Base", "Fighting Force: U.S. Army Brigade", "Tripwire: 240 U.S. Soldiers", "Offshore: Stockpiling Supplies", "Offshore: U.S Navy Destroyer", "Transient Demonstration: Bomber/Fighter Patrols"),
  mean=c(3.958, 3.844, 2.109, 3.213, 3.2, 2.978),
  se=c(0.114, 0.121, 0.13, 0.152, 0.135, 0.126)) 

tiff("Figure4_Capability.tif", res=100, height=10, width=10, units="in")
CapabilityData %>%
  arrange(mean) %>%
  mutate(name = factor(name, levels=c("Fighting Force: Permanent U.S. Base", "Fighting Force: U.S. Army Brigade", "Offshore: Stockpiling Supplies", "Offshore: U.S Navy Destroyer", "Transient Demonstration: Bomber/Fighter Patrols", "Tripwire: 240 U.S. Soldiers"))) %>%
  ggplot(aes(x = name, y = mean)) + geom_point() +
  geom_errorbar(aes(ymin = mean - 1.96 *se, ymax = mean + 1.96 * se), width = .2) +
  theme_minimal(base_size=12, base_family="EB Garamond") + xlab("Reassurance Measure") + ylab("Mean Perceived Capability") +
  theme(axis.text.x = element_text(hjust = 1)) + coord_flip()
dev.off()


## Figure 5

CapabilityData2 <- data.frame(
  Measure=c("Fighting Force: Permanent U.S. Base", "Fighting Force: U.S. Army Brigade", "Tripwire: 240 U.S. Soldiers", "Offshore: Stockpiling Supplies", "Offshore: U.S Navy Destroyer", "Transient Demonstration: Bomber/Fighter Patrols"),
  CapabilityMean=c(3.958, 3.844, 2.109, 3.213, 3.2, 2.978),
  CapabilitySE=c(0.114, 0.121, 0.13, 0.152, 0.135, 0.126)) 

LikelihoodData2 <- data.frame(
  Measure=c("Fighting Force: Permanent U.S. Base", "Fighting Force: U.S. Army Brigade", "Tripwire: 240 U.S. Soldiers", "Offshore: Stockpiling Supplies", "Offshore: U.S Navy Destroyer", "Transient Demonstration: Bomber/Fighter Patrols"),
  LikelihoodMean=c(3.978, 3.93, 3.395, 3.311, 3.205, 3.209),
  LikelihoodSE=c(0.165, 0.145, 0.129, 0.137, 0.14, 0.149)) 

z <- merge(LikelihoodData2, CapabilityData2, by="Measure")

tiff("Figure5.tif", res=100, height=10, width=10, units="in")
ggplot(z, aes(x=CapabilityMean, y=LikelihoodMean, group=Measure)) + geom_point(aes(shape=Measure)) + 
  geom_errorbar(aes(ymin = LikelihoodMean - 1.96 *LikelihoodSE, ymax = LikelihoodMean + 1.96 * LikelihoodSE), width = .05) +
  geom_errorbar(aes(xmin = CapabilityMean - 1.96 *CapabilitySE, xmax = CapabilityMean + 1.96 * CapabilitySE), width = .05) +
  theme_minimal(base_family="EB Garamond") + xlab("Mean Perceived Capability") + ylab("Mean Perceived Likelihood that U.S. Will Aid Ally in Defense") +
  theme(axis.text.x = element_text(hjust = 1), legend.position = "bottom") 
dev.off()
   


##APPENDIX##

##Demographics

##Country
CountryTable <-table(data$Country)
prop.table(CountryTable)

##Government Experience
GovExperienceTable <-table(data$GovExperience)
prop.table(GovExperienceTable)

##Current Affiliation
AffiliationTable <-table(data$CurrentAffiliation)
prop.table(AffiliationTable)

##Gender 
GenderTable <-table(data$Gender)
prop.table(GenderTable)

mean(data$Age)




    
#Correlation Matrix: Capability vs. Resolve vs. Safety
data$capability <- (data$Base_Capability + data$Brigade_Capability + data$Company_Capability + data$Stockpile_Capability + data$Destroyers_Capability + data$Flybys_Capability) / 6
data$resolve <- (data$Base_Resolve + data$Brigade_Resolve + data$Company_Resolve + data$Stockpile_Resolve + data$Destroyers_Resolve + data$Flybys_Resolve) / 6
data$safety <- (data$Base_Safe + data$Brigade_Safe + data$Company_Safe + data$Stockpile_Safe + data$Destroyers_Safe + data$Flybys_Safe) / 6

df1 <- as.data.frame(cbind(data$capability, data$resolve, data$safety))
df1$Capability <- df1$V1
df1$Resolve <- df1$V2
df1$Safety <- df1$V3
df1$V1 <- NULL
df1$V2 <- NULL
df1$V3 <- NULL

library(corrplot)
source("http://www.sthda.com/upload/rquery_cormat.r")
rquery.cormat(df1)

dev.copy(png,'CorrelationMatrix.png')
dev.off()
